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Abstract. We study the zero-range process on a simple diamond lattice with open 
boundary conditions and determine the conditions for the existence of loops in the 
mean current. We also perform a large deviation analysis for fluctuations of partial 
and total currents and check the validity of the Gallavotti-Cohen fluctuation relation 
for these quantities. In this context, we show that the fluctuation relation is not 
satisfied for partial currents between sites even if it is satisfied for the total current 
flowing between the boundaries. 
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1. Introduction 

In the last decades much effort has been put into understanding and modelling non- 
equilibrium systems, which find applications in fields ranging from biology to finance 
[H, [2]. Among the various models which have been proposed to study such real-life 
processes, stochastic interacting particle systems (interacting Markov systems) have en- 
joyed particular success [3j. In this class, the zero-range process (ZRP) is a well studied 
lattice gas model offering many applications and the possibility of obtaining analytical 
results. Introduced in 1970 [1], one of the reasons the ZRP gained interest was because 
it can show a phase transition from a fluid to a condensed state |5] . The ZRP has been 
extensively studied with both periodic and open boundary conditions in one dimension 
P,[7], and some variants involving junction topologies have also been introduced, see e.g. 
[SI [9]. Furthermore, currents in a closely related model have recently been considered 
on more general networks ^Q\, which may give some insight into expected effects for 
manmade networks such as traffic on roads or the internet. 



For extended ZRPs defined on two-dimensional or three-dimensional lattices, parti- 
cle currents can flow in principle in loops within the bulk of the system. To illustrate this 
possibility, we study a simple variant of the ZRP defined on a diamond-shaped lattice 
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with open boundaries. For this model, we are able to obtain analytical expressions for 
the mean current flowing between specific sites of the lattice which show that the model 
has two regimes of mean particle currents: a unidirectional regime where particles flow 
in the same direction through different lattice branches and a loop-current regime where 
a mean current flows around the diamond. 

In addition to mean currents, we also study fluctuations by calculating the proba- 
bility distribution of particle currents and its associated large deviation rate function, 
which plays a role similar to thermodynamic potentials in equilibrium systems [TTl [12] . 
Recently much experimental and theoretical attention has been devoted to the study 
of certain fluctuation symmetries [131 El; which may also be observed experimentally, 
e.g. [HI m]. Our simple model allows us to study explicitly the joint probability of 
observing a given current on different lattice branches. We obtain analytically the large 
deviation functions for currents across different bonds and hence gain understanding of 
the role of current loops. These results for current loops are expected to be important 
in testing fluctuation symmetries for higher dimensional systems, such as those reported 
recently by Hurtado et al. [T7] . 

The remainder of the paper is structured as follows. In section [21 we define the ZRP 
on the diamond lattice and calculate its stationary state. In section [3l we discuss the 
appearance of unidirectional and loop-mean-current regimes in the system. In section 
[H we calculate the joint particle current fiuctuations of the upper and lower branches of 
the lattice and test the well known Gallavotti-Cohen fiuctuation relation for partial and 
total currents. Finally, in section [5l we discuss our results and their possible implications 
for other models. 

2. Zero-range process on a diamond lattice 

2.1. Definition of model 

The ZRP is a lattice-based many-particle model in which, as the name suggests, par- 
ticles interact only with other particles at the same site. Additionally, particles are 
allowed to accumulate to any non-negative number on each site of the lattice. Here, we 
study the ZRP with open boundary conditions on a diamond lattice as shown in figure [H 

The dynamics of the particles on this lattice is defined in continuous time such 
that the topmost particle on each site hops to an adjacent site after an exponentially 
distributed waiting time. More precisely, particles hop clockwise around the diamond 
with rate pWn and anti-clockwise with rate qWn, where the interaction between particles 
on each site is taken into account by the term Wn which depends exclusively on the 
occupation number n of the departure site. Given the symmetry of the system, without 
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Figure 1. Diamond array. The z^'s represent the fugacities of sites 1-4. The quantities 
a, 7, S, j3, p and q are the hopping rates. 



loss of generality we will assume from now on that 
p > q. 



(1) 



Moreover, we allow particles to enter and leave the boundary sites with probability rates 
a and (3wn for site 1, S and ■yWn for site 3. 

Based on the one-dimensional open boundary ZRP studied in [7j we expect the sys- 
tem to be driven out of equilibrium for a/7 7^ S/P but that some choices of parameters 
may lead to a boundary condensation phenomenon in which particles accumulate on 
one of the sites 1 or 3. We shall return to this point later. 



In the quantum Hamiltonian formalism [18], one defines a probability vector 

\P) = J2Pin)\n) (2) 

n 

where \n) is a basis vector for the particle configuration n = (ni, 71-2, ""-3, and P{n) 
the probability of that configuration. Then the time evolution of the system is described 
by the master equation 
d\P) 



dt 



-H\P). 



(3) 



Here the matrix if, or Hamiltonian, is the stochastic generator of the system. To 
explicitly write the Hamiltonian of our system we define the creation and annihilation 
operators on site i by 

\ / Mil 





/ 
1 
1 



V 



and a,- 









W2 





(4) 
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respectively. With the additional definition of the diagonal matrix di = wjdj^k, the 
Hamiltonian of the model shown in figure [1] is written 

- H = a {a^ -l) +-f {a^ - di) + 6 {a^ - l) + /3 (a^ - 4) 

4 

+ J2p -dk)+q K«fe+i - (5) 

fc=i 

where the index sums are taken modulo 4. 
2.2. Steady state 

We are interested in finding the non-equilibrium stationary state \P*) of our system. 
By definition, this probability does not change in time and is therefore, such that 

H\P*) = 0. (6) 

It has been shown that for the ZRP on any lattice geometry the steady state factorises 
as the tensor product (see, e.g., [THj) 

\P*) = |P*(ni)) ® \P;in2)) ® ... ® |P;K)). (7) 

Here each site's probability distribution vector is given by 

\P:M) = J2Pk{nk)\n,), (8) 

with Pkiuk) the probability of finding rik particles on the /cth-site and \nk) the 
corresponding configuration vector. Furthermore, Pk{nk) is given in terms of the site's 
fugacity Zk and the interaction term w„ by 

Pk{nk) = Z~'zl^l[w;\ (9) 

i=l 

where Z is the grand canonical partition function 

oo j 
j=0 1=1 

It now becomes obvious that the choice of the interaction Wn will be refiected in the 
existence of the partition function. For instance, if we let 

lim Wn = (11) 

n— >oo 

with K a constant, then we have to make sure that z < n in order for the infinite sum 
of equation f ITU]) to converge and the partition function to exist. Hopping parameters 
leading to z > k, correspond physically to a growing condensate. In the remainder of 
this paper we will assume, unless stated otherwise, an unbounded interaction rate Wn 
for the particles, i.e., 

lim Wn = oo. (12) 



n—>-oo 
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This guarantees that the system has a non-equihbrium stationary state without con- 
densation. 



It can be shown that the creation and annihilation operators act on the stationary 
state eigenvector of the Hamiltonian as a^\Pk) = z'^'^dk\P^) and ak\P^) = Zk\P^) 
respectively. Here we have made explicit that the operators corresponding to the fcth- 
site act only on the probability distribution vector of the same site. Thus, applying the 
Hamiltonian (j5]) to the stationary state vector \P*) yields the expression 

H\P*) = - [{a - + qz2 - pzi + pz4 - qzi) z^^di 
+ {pzi - qz2 + qzs - PZ2) ^^^2 
+ (5 - l^z^i + pz2 - qz3 + qZi - pz^) z^^d^ 
+ {qzi - pz4 + pzs - qz4) z^^d^ 

-{a-^zi + 6-Pz3)]\P*) (13) 

and for \P*) to be the required stationary state the coefficients of the matrices di must 
vanish. Solving the resulting system of equations leads to the fugacities 

(p + q)a(3 + {p"^ + q'^){a + 5) 



Z2 

zz 



(p + g)/37 + (p2 + g2)(/3 + 7)' 
pa/3 + q^b + (p^ + + b) 

(p + g)/37 + (p2 + g2)(/3 + 7)' 
(p + g)/37 + (p=^ + g2)(/3 + 7)' 



Note that with these solutions, one verifies that a — 721 + 5-/3^3 = 0, which is consistent 
with the stationary state being the eigenvector with eigenvalue zero. 

With these fugacities, we can calculate the mean time-averaged current ja,h which 
measures the average number of particles jumping from site a to site h per unit time. 
This current is expressed in terms of the fugacities as 

jfc.fe+i = VZk — qzk+i = —jk+i,k- (15) 

Here once again, k G {1,2,3,4} and modulo 4 applies over the subindex addition. 
Similarly, for the boundary sites the mean current is 

= a - 721, 

jn = Pzs - 5, (16) 

with the convention of a positive fiow direction from the left to the right boundary of 
the lattice. 
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Since we have assumed the system is in a stationary state without condensation, 
the mean currents must satisfy 

3L = jl,2 + jl,4 = J2,3 + J4,3 = JR (17) 

by particle conservation. From this and the assumption that p > q, we then find that 
the mean current at the left boundary is positive if 

a-7^i>0, (18) 

which is equivalent to 

ip' + q') (a/3-75) 



(p + g) /37 + (P' + g') (/3 + 7) 



> (19) 



or 



a/3 - 75 > 0. (20) 

This condition implies a net current left-to-right through the system. If conditions ([T]) 
and fl20l) are satisfied, it follows that the largest fugacity is Zi. So returning to the 
discussion of a bounded interaction Wn, the consistency condition for existence of a 
stationary state without condensation is Zi < k. 

3. Particle current loops 

An interesting feature of the model we are studying is that it shows a transition from a 
unidirectional-mean-current to a loop-mean-current regime as the parameters are var- 
ied - see figure [2l An intuitive way to see this change, is to consider specific parameter 
values. On one hand, when we have symmetric hopping rates p = g = 1, the mean 
current fiows clockwise through the upper branch and ant i- clockwise through the lower 
branch, i.e., the unidirectional regime. On the other hand, when we have completely 
asymmetric hopping rates p = 1 and g = 0, the mean current fiow is forced to go around 
the diamond, i.e., the loop-current regime. Since the mean current through the lower 
branch changes from anti-clockwise to clockwise direction as the regime changes from 
unidirectional to a loop, one way to determine where the transition occurs is to calculate 
the parameters at which the mean current ji^4 between sites 1 and 4 vanishes. 
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As a first step to characterize this change of regime in more detail we now consider 
the special case of 7 = 5 = 0, which corresponds to having only injection of particles 
at the left-end of the lattice and only depletion at the right-end. This allows us to 
understand the change of regime with fewer parameters and get some intuition for the 
cases with non-zero hopping rates 7 and 6. 



From equations (|T4|) and (|T5|l the condition for a vanishing current ji^4 is satisfied 



when 



B 



(21) 



where we have defined the ratios Q = q/p and B = P/p. The curve defined by the 
equation above gives us the exact location of the regime change from positive to negative 
average current as shown in figure [31 As expected, if Q = (i.e., g = 0) it does not 
matter how large the extraction rate B (or /3) is, the system is always in the loop 
regime, whereas if Q = 1 (i.e., for symmetric hopping rates) the system is always in 
the unidirectional regime. Notice also that equation ( 12T|) has no dependence on the 
injection parameter a which means, that for < Q < 1 we can control completely the 
system regime just by changing the extraction rate /3. For large /3, the system favours 
the unidirectional regime whereas, for low /3, a smaller fraction of the particles on site 
3 can leave the system and a loop current is therefore more likely. 




Figure 3. Regimes of the average current ji.4. Red region ji 4 < 0. Blue region 
jl,4 > 



In the general case where we admit injection and extraction of particles from both 
boundary sites, we can still compute the mean current ji^4: 

-■ ^ (g + 6) {p^ + q^) {q-p) + al3q^ - -jSp^ , . 

{p^ + q^W + l) + {p + q)h ■ ^ ^ 

We can check again that for the special cases of symmetric hopping rates (i.e., p = q = ^) 
the current is positive and the system is in the unidirectional regime, whereas for totally 
asymmetric rates (i.e., p = 1 and q = ^) the current is negative and thus, in the 
loop-current regime as expected. Moreover, for fixed p > g > 0, we can see that the 
mean-current regime can be chosen by changing the boundary rates where, just as in 
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the special case discussed above, large /3 favours the unidirectional regime and small (3 
the loop-current regime. Similarly, for fixed boundary conditions satisfying a/3 > 7^, 
we can choose the mean-current regime by varying the bulk hopping rates p and q. A 
further remark about the existence of these two regimes, is that if we were to consider Wn 
bounded we know that for some choices of the hopping rates the system would undergo 
condensation (and have no steady state) meaning, for example, that not all of the phase 
plane in figure [3] would be accessible. 

4. Particle current fluctuations 
4.1. Large deviations 

We now complement our results for the mean currents by studying their fluctuations. 
To be specific, we are interested in calculating the probability distributions p{ja,b,t) of 
the time-averaged currents ja,b between sites a and b in the lattice. In the long-time 
limit, we expect this distribution to follow the large deviation principle 

Pija,b, t) ~ exp {-te{ja,b)) , (23) 

with rate function (RF) e{ja,b)- Here "~" denotes asymptotic equality in the limit of 
large time. One way to obtain the RF is to calculate the log-scaled cumulant generating 
function (SCGF) 

e(A) = lim --log(e-^*^«-''). (24) 

t— >oo t 

Indeed, if the SCGF is continuous and differentiable, then it is known that the RF is 
obtained from the SCGF via the Legendre transform [11] as 

e{ja,b) = max {e(A) - \ja,b} ■ (25) 

A 

To calculate the SCGF we need to modify the Hamiltonian ([5]) to count particle 
jumps. This is done by multiplying the terms accounting for the transfer of particles 
from sites a to 6 (and vice versa) by the exponential factor e^^, see for example [13] . 
Then for instance, to measure the fluctuations of ji,4 the required modified Hamiltonian 
is: 

- H = a{at -1) -^'yla^ -di) + 6 (a^ - l) + (3 {a^ - 4) 

4 

+ - 4) + g {at^+.e-''"-^ - 4+1) . (26) 

k=l 

The SCGF can be identified with the lowest eigenvalue of the Hamiltonian fl26l) . which 
in general is different from zero. Then, we find the desired RF of the current ji.4 via 
the Legendre transform fl25l) . 

4-2. Joint probability of current fluctuations 

Correlations between particle currents flowing through different bonds can be studied 
by considering a two-parameter SCGF. For the diamond lattice, it is interesting for 
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example to study the current flowing between sites 1 and 2 and the current flowing 
between sites 1 and 4 simultaneously. The joint SCGF that characterize the joint 
distribution p (ji,2, ji,4, ^) of these currents is calculated analytically by modifying the 
system's Hamiltonian similarly to described above. Particle jumps between sites 1 and 
2 are counted with parameter Ai^2 and between sites 1 and 4 with parameter Ai^4. This 
calculation leads to the two-parameter SCGF, 

e[ 1,2, l,4j - Q4^2>^i,2_^_e^>^iA-e^i,2+HA(l+Q4+B{Q+l){Q^+l)+G{Q+l)(Q^+BQ+B+l)) 

(a+<5)(Q4e2^1,2-e^l,2+^l,4(Q4+i^+g2Ai_4^ 
+ Q4e2Ai,2+e2Ai^4_gAi,2+Al,4(l+Q4+B(Q+l){Q2 + l)+G(Q+l)(Q2+BQ+B+l)) ^ ^ 

where we have defined Q = q/p, B = P /p and G = '~f/p. For the special case 7 = 5 = 
the SCGF reduces to 

p(\ \ I B(l+Q)(e^^-^+e'i.2Q2) \ , . 

e{Al,2, Al,4j - a + ^2Ai,4+e2Ai,2Q4_eAi,2+Ai,4(Q4 + i + (i+Q)(i+Q2)5) ) ' l^^J 

The RF e(ji^25 ji,4) associated with the joint distribution p {ji,2, ji,'i,t) is obtained from 
this SCGF via the double Legendre transform: 



e(ji,2, ji,4) = max {e(ji,2, ji,4) - ji,2Ai,2 - ji,4Ai,4} • (29) 

-'ll,2,Ai,4 

The numerical evaluation of this transform is shown in figure |H 




Figure 4. Joint probability rate function with parameters a — 1/2, G = 0, S = 0, 
B = I and Q = 1 (i.e., unidirectional regime). 

The RF of figure H] is not defined for all currents, since for 7 = 5 = it follows that 
Ji,2 + Ji,4 > 0. For the parameters used in figure HJ we also know from our analysis of 
section [3] that each of the mean currents ji^2 and ji^4 corresponding to the minima of 
e(ji,2, Ji,4), are strictly positive, which means that we are in the unidirectional regime. 
What the joint RF shows is that, despite being in this regime, it is possible to have 
loop current fluctuations. In particular, for a fixed large positive current ji,2, the most 
probable value of ji,4 is negative (as this minimises the RF), implying that loop current 
fluctuations are more likely than unidirectional current fluctuations in this case. This 
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behaviour is also seen for other values of parameters. 

Note that from the joint SCGF we can quickly recover the SCGF of either of the 
two possible partial currents or the total current simply by making the correct selection 
of parameters Ai^2 and Ai,4. We shall do this in the following subsections. 



4.3. Partial current fluctuations 

A possible use of the two-parameter SCGF obtained above is to analyse the partial 
currents ji,2 or ji^4. Since the particle flow from site 1 to site 4 allows us to observe 
both positive and negative mean currents (see section 3), we choose to analyse this bond 
more closely. To do this we set the parameters Ai,2 = and Ai,4 = A and the resulting 
function is 

We can further simplify this function by considering the case 7 = 5 = in which 
we already studied the change of regime from unidirectional to loop- mean-current. The 
resulting SCGF and the numerical RF are shown in figure [51 Here we can observe again 
how the choice of the parameter B affects the mean current flow through the lattice. 




Figure 5. SCGF (above) and RF (below) for current ji^^ with parameters a — 1/2, 
Q = 1/2 and S = 3/2 (red), (black), B 9/2 (blue). 

We can see directly from the SCGF whether the choice of parameters corresponds 
to a unidirectional or loop regime by looking at the slope of the function at A = 0. This 
is because, from large deviation theory, we know that the first derivative of the SCGF 
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at zero corresponds to the mean current between the sites we are looking at. Hence, 
a positive slope implies a positive current and the unidirectional regime, whereas a 
negative slope implies the loop regime. At the level of the RF, this translates into having 
a positive or negative minimum, respectively, which determines the mean current. From 
the full form of the RF, what can be seen again is that fluctuations of currents having 
a sign opposite to the mean current have a non-zero probability to be observed. This 
means concretely that loop current fluctuations can be seen in the unidirectional-current 
regime and vice versa. 



4.4- Total current fluctuations 

Another possibility, is not to differentiate between the two lattice branches but simply 
to measure the total flux of particles through the cross-section between site 1 and sites 
2 and 4, i.e., j = ji,2 + ji,4- Once again, the joint SCGF (!27j) reduces to a single variable 
function but now by taking Ai,2 = Ai_4 = A: 

n u (1 - c-') iQ' + 1) [aB - SGe') 
<^>-<^>^)= BiQ' + l) + GiQ' + BQ + B + iy 
The RF e(j) associated with the probability distribution of j can be obtained 
analytically for this SCGF, resulting in 

-(p^+g^)M+7.5) 



e(j) 



(p2+'?2)(/3+7)+(P+f/)/37 
I a/ (j(p^+g^)/3+i(p^+9H(p+g)/3)7)V4(p2+g2-)2^^ 

(pH9")(/3+7)+(P+9)/37 
i(p/37+p2(/3+^)+g(q/3+(g+/3)^)) 

j(p/37+p^(/3+7)+g(g/3+(g+^)7)) 
2(p2+g2)7<5 



X log 



"t" 2(p2+g2)7<5 )• y'^'^) 

As expected, for the case of unbounded Wn the same RF is found when measuring 
current fluctuations of j2,3 + j4,3 or either of the boundary currents Jl or Jr. We shall 
use equation ( 132|) in the next subsection to show that fluctuations of total currents 
possess a particular symmetry property that fluctuations of partial currents do not 
have. 



4-5. Fluctuation symmetries 

In the previous subsections we have seen that the large deviation RF provides rich in- 
formation about stochastic models beyond the level of mean values. The study of such 
fluctuations allows us to find symmetries in the dynamics of a system rather than just 
spatial symmetries. In particular, the Gallavotti-Cohen fluctuation relation (GCFR) 
provides a well-known symmetry for the entropy production rate, which can be trans- 
lated into a symmetry for the current RFs in interacting particle systems [131 ISHl HH |22] . 
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The GCFR in the context of particle current fluctuations can be written as 



-Ejt 



. .,^ ^ (33) 

and makes a connection between the long-time probability of observing a time-averaged 
current j and the probability of observing a current —j of the same magnitude in 
the reversed direction. Here, E is a. constant which is commonly interpreted as an 
equilibrium-restoring field. Since we are assuming that currents follow a large deviation 
principle as in equation (123|1 . relation (133|1 can also be written in terms of only the RFs 
as 

- e(-j) = -Ej. (34) 

Now that we know the probability rate function of observing particle currents 
measured through different bonds, we would like to investigate if our rate functions 
satisfy the GCFR. In our model it can be straightforwardly shown analytically that the 
RF (132|) for the fluctuations of the total current, j = ji,2 + ji,4, obeys the GCFR (l34l) 
with 

'a/3' 



E = log 



7(5 



(35) 



Interestingly, due to the cyclic arrangement of the bulk hopping rates in this model, the 
equilibrium-restoring field E turns out to be independent of the rates p and q and equal 
to the expression for a single-site ZRP with open boundaries 



On the other hand, for partial current fluctuations, the GCFR is not satisfied. 
This can be clearly seen in figure [6] where we plot e(ji,4) — e(— ji^4) for the numerical 
results obtained in section 4.3 and observe that is not linear in ji 4. However, we can 
see that for small ji^4 the relation is approximately satisfied. 




Figure 6. GCFR breakdown for partial current Ttie blue curve corresponds 

to the numerical values of e(— ^1^4) — e{ji^i) with parameters Q = 1/2, a ~ 1/2, 
G = 5 = and B = 9/2 (unidirectional regime) whereas the straight line corresponds 
to e{—jiA) — e{jiA) = EjiA with field E = 0.563. 

Our observations are consistent with recent results for other two dimensional 
models, both classical and quantum, which highlight the breakdown of the GCFR for 
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partial currents [23] 1^. We emphasize here that the crucial point is not the existence 
of a loop-mean-current regime but that for some partial current fluctuations away from 
the mean, the most likely realisation involves loop current flow. This is confirmed by 
studying a variant of our model in which the rates on the lower branch are reversed (see 
figure [7j) and thus the possibility of a loop-mean-current is eliminated. Performing a 
similar analysis as for the original diamond model one again observes loop fluctuations 
in the joint-probability RF and finds that the GCFR is satisfied for the total current, 
this time with E = but not for partial currents. 




Figure 7. Model variant with reversed hopping rates on lower branch 
5. Conclusion 

We have performed a current fluctuation analysis of the zero-range process on a dia- 
mond lattice with open boundaries. For rates which impose a preferred direction around 
the diamond, we demonstrated the possibility of two different mean current regimes de- 
pending on the boundary parameters and the asymmetry of the bulk hopping rates. 
Moreover, for the case in which particles are only injected on one side of the diamond 
and only removed from the other we proved that, regardless of the injection rate, we 
could control the regime of the system just by varying the extraction rate. 

Analysing the current fluctuations via the two-parameter scaled cumulant generat- 
ing function, we could study the joint probability distribution of the partial and total 
currents flowing between sites 1 and 2 and between sites 1 and 4 of the lattice. Sig- 
nificantly, we saw from the joint current rate function that both unidirectional- and 
loop-current fluctuations may be observed whether the mean current is unidirectional 
or in a loop. 

From our analysis we also confirmed that to observe the Gallavotti-Cohen fluctua- 
tion relation we need to measure total currents. This point should obviously be taken 
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into account when testing fluctuation symmetries in higher dimensional systems in ex- 
periment or simulation. Indeed, the recent isometric fluctuation relation [17] which is a 
generalised symmetry for higher dimensions, is also concerned with global rather than 
local currents. 



Possible extensions of this work would be to consider larger higher dimensional 
lattices or more complex geometries. For example, we could characterize a chain of 
diamonds similar to the one studied in [26] for the Ising model. We emphasize that our 
results hold only for the case of an unbounded site interaction w„; it would be interesting 
to extend the analysis to the case of bounded Wn where one expects the formation of 
"instantaneous condensates" and a breakdown of the fluctuation symmetry even for 
total currents [23l |27l EH]- In all of these models we expect loop-current or vortex 
fluctuations to play an important role. 
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